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SECTION  1 


INTRODUCTION 

During  the  period  of  this  contract  Berkeley  Research  Associates 
(BRA)  was  tasked  to  address  the  problem  of  structure  creation  and  con¬ 
vection  In  the  Intermediate  and  low  altitude  regime.  More  specifically 
the  effort  was  directed  toward  finding  alternative  methods  by  which  the 
nonlinear  dynamics  of  structure  formation  and  convection  could  be  studied. 
If  such  approaches  proved  adequate,  then  the  Information  generated  could 
be  used  as  models  In  various  systems  codes  where  the  knowledge  of  both 
neutral  fluid  structure  and  plasma  structure  Is  very  Important  to  evalua¬ 
tion  of  various  systems. 

The  effort  took  two  directions.  One  was  to  Investigate  the  applic¬ 
ability  of  current  microstructure  algorithms  to  the  Intermediate  altitude 
structure  and  convection  problem.  The  second  direction  was  the  construc¬ 
tion  and  testing  of  a  vortex- In-cell  approach  to  investigate  this  technique 
ability  to  handle  vortex  formation  and  torusing  dynamics  in  multiburst 
situations.  Further  these  techniques  offered  the  possibility  of  studying 
the  production  of  very  small  scale  structure  evolution. 

During  the  course  of  this  contract  a  2-D  vortex- in-cell  code,  VINCE, 
was  constructed  and  tested.  Methods  by  which  VINCE  could  be  generalized 
to  3-D  were  considered  and  found  to  be  something  we  could  accomplish. 

The  high  altitude  microstructure  algorithm  was  considered  as  a  model 
for  a  possible  mid-altitude  structure  algorithm  and  avenues  of  further 
research  were  laid  out. 

In  the  following  sections  a  discussion  of  the  results  obtained  in 
this  third  of  a  man  year  effort  will  be  presented.  In  Section  2  a  pre¬ 
sentation  will  be  made  of  the  numerical  results  to  date.  These  results 
Include  some  of  the  possible  nonlinear  dynamics  that  might  occur  in  the 
multiburst  environment.  In  Section  3  a  discussion  of  some  concept  con¬ 
cerning  a  possible  structure  convection  model  will  be  presented.  Section 
4  will  contain  a  brief  discussion  of  the  simulations  accomplished  during 
this  contract. 
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SECTION  2 
SIMULATION  METHOD 


As  mentioned  previously,  one  of  the  tasks  was  to  address  the  problem 
of  structure  creation  In  the  low  to  Intermediate  altitude  regime.  In 
this  context  Intermediate  means  bursts  of  up  to  90  km  or  so  and  low  means, 
from  our  standpoint,  altitudes  above  which  ground  effects  need  to  be 
considered.  One  of  the  most  obvious  mechanisms  by  which  structure  is 
created  Is  the  toruslng  of  the  fireball  which  Is  a  form  of  the  Rayleigh- 
Taylor  Instability.  Less  obvious  but  maybe  as  Important  Is  the  further 
evolution  of  these  torused  regions  when  they  are  In  the  flow  field  created 
by  other  bursts,  or  the  convection  of  debris  material  in  the  flow  fields 
of  this  late  time  fireball.  Here  late  time  implies  after  the  original 
shocks  have  left  the  immediate  vicinity  of  the  problem. 

A  tool  that  might  provide  Information  of  the  sort  just  mentioned 

Is  a  vortex- in-cell  code.  The  basic  algorithm  was  first  explained  and 

used  by  Meng  and  Thomson.1  It  makes  use  of  the  fact  that  physical  problems 

such  as  shear  flow,  flow  with  a  sharp  Interface  separating  two  fluids 

of  different  densities  or  regions  of  different  viscosities  can  be  handled 

by  noting  that  generation  and  convection  of  vortlcity  is  the  primary 

quantity  of  interest.  For  example,  because  of  the  density  difference 

across  the  Interface,  vorticity  is  generated  along  it  by  the  interaction 

of  the  pressure  gradient  and  density  gradient.  The  resulting  motion 

consists  of  essentially  two  processes:  the  creation  of  a  vortex  sheet 

and  the  subsequent  mutual  Induction  of  different  portions  of  the  sheet. 

The  method  treats  vorticity  as  a  basic  quantity  of  fluid  motion  to  be 

handled  in  a  code  much  the  same  way  charged  particles  are  treated  in 

a  plasma  particle  simulation.  The  algorithm  as  discussed  by  Meng  and 

Thomson1  simply  follows  the  interface  between  the  hot  fluid  and  the  cold 

fluid  background,  computes  the  development  and  convection  of  the  vorticies 

which  in  turn  delineate  the  dynamics  of  the  Interface.  These  calcula- 
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tlons  can  be  done  in  both  viscid  and  Inviscid  situations.  ’  ’  As  will 
be  shown  later,  the  code  calculates  the  torusing  of  a  rising  fireball 


or  bubble.  This  formulation  can  also  address  smaller  scale  structure  such 
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as  the  Kelv in-Helmholtz  instability  and  the  Saffman-Taylor  ’  instability, 
both  of  which  may  play  a  role  in  the  structuring  of  debris  material  and 
perhaps  the  rising  fireball. 

The  algorithm  Itself  uses  a  very  simplified  version  of  the  Navier- 
Stokes  equation.  Only  enough  of  the  nonlinearity  is  retained  to  provide 
the  structure  evolution  of  Interest.  One  begins  with  the  momentum  equation 
for  an  incompressible  fluid 

*  -vp/p  +  n  v2u  +  g  (1 

where  u  is  the  velocity,  p  is  the  pressure,  p  the  fluid  density, 
n  the  kinematic  viscosity  and  g  is  the  acceleration  of  gravity. 

Vortlcity  is  defined  as 

l  *  v  x  u  (2 

and  Equation  (1)  becomes 

§  »  v0)  x  vp  +  n?2c  •  (3 

The  algorithm  can  be  set  up  to  do  a  viscid  calculation  or  this  equation 
can  be  further  reduced  by  making  the  inviscid  assumption.  One  then  has 
the  operative  equation  used  in  the  simulations  to  be  shown  in  this  report. 
Equation  (3)  becomes 

.  -V  (!)  x  7p.  (4 

A  further  assumption  can  be  made  in  the  case  of  torusing  bubbles  that 
still  retains  the  nonlinearity  of  the  flow  field.  This  is  the  assumption 
that  the  pressure  can  be  treated  to  only  first  order  then 
vp  *  P09  +  0(Vp/pQ).  Equation  4  then  can  be  written  as 
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Here  we  are  working  in  the  x-z  grid  where  g  *  -g  n  .  As  expected  the 


vortlcity  is  directed  out  of  the  plane.  From  Equation  (5)  one  sees  that 


the  vorticity  is  generated  only  on  the  interface  between  the  two  fluids. 
The  total  circulation  for  an  ith  element  of  the  fluid  is  given  by 


dx'dz'. 


The  evolution  of  the  ith  circulation  element  then  becomes 


(6) 
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where  pQ  =  -|(p+  +  p_)  and  p+  and  p_  are  the  densities  across  the 
interface  and  az^  is  the  vertical  dimension  of  the  fluid  element.  This 
formalism  leads  naturally  to  the  concept  of  tracking  various  fluid  elements 
to  get  the  total  dynamics  of  the  system.  The  standard  method  of  tracking 
this  element  is  to  consider  the  effects  of  all  other  vortex  elements  on  the 
one  being  considered.  This  is  the  Green's  function  formalism.  An  alternate 
method  is  to  express  the  velocity  as  a  stream  function.  Meng  and  Thomson1 
compared  both  methods  and  found  the  latter  formalism  much  faster  hut  not 
quite  as  accurate  on  very  small  scale  phenomena.  Because  of  our  interest 
in  a  possibly  large  number  of  cases  the  later  approach  was  attempted. 

If  one  assumes  that 


U  *  V  X  £  (8) 

where  v  •  iji  =  0  (an  automatic  requirement  in  2  dimensions),  then  the 
stream  function  satisfies  the  Poisson  equation 

v2$  *  -5  (9) 

where  fl  i  £  5.  S(x  -  x. ),  a  summation  to  an  Eulerian  grid.  The  system 
i  1  9  1 

is  now  closed  because  one  knows  the  new  stream  function  and  hence  the  new 

velocities  and  positions  from  Equations  8  and  9.  Equation  7  then  advances 

the  circulation  and  the  process  proceeds.  The  algorithm  is  written  very 

much  like  an  electrostatic  particle  code  with  a  time-centered  leap  frog 

algorithm,  particle  interpolation  and  weighting.  Similar  formalisms  have 

12  3 

been  developed  for  the  viscid  calculation.  ’  * 


At  this  juncture  we  have  a  vortex-in-cell  formalism  where  we  are 
solving  an  Incompressible  fluid  problem  with  the  assumption  of  a  linearized 
pressure.  What  can  such  a  simple  model  produce?  One  of  the  most  obvious 
and  Interesting  nonlinear  result  is  the  Raylelgh-Taylor  instability  which 
in  the  case  of  a  buoyantly  rising  bubble  is  the  torusing  evolution.  We 
have  performed  a  variety  of  such  simulations  to  test  out  the  code  and 
also  get  a  glimpse  of  the  possible  applications.  It  should  be  recalled 
that  this  algorithm  Is  not  intended  to  replace  detailed  hydro  codes. 

However,  it  will  allow  the  Investigation  of  the  nonlinear  interaction 
of  flow  fields  such  as  one  might  expect  in  the  late  time  evolution  of  a 
multiburst  situation. 

One  of  the  first  calculations  that  we  performed  with  the  new  code 
called  VINCE  was  to  calculate  the  evolution  of  a  bubble  into  a  torus  or 
vortex  pair.  As  mentioned  by  Meng  and  Thomson1,  one  can  get  quantitative 
estimates  of  the  rise  rate  and  separation  of  the  pair  by  observing  as 

4 

Turner  did  that  the  circulation  of  each  vortex  approaches  a  constant 
value  after  vortex  or  torus  formation.  After  the  vortex  formation  the 
density  along  a  path  threading  the  center  of  the  vortex  is  essential lly 
constant  and  equal  to  the  ambient  value  and  dr/dt  +  0.  For  a  constant 
r  the  rise  velocity  varies  inversely  as  the  separation  of  the  vortex 
pair;  u  «  1/R.  The  upward  momentum  increases  at  a  constant  rate 
d(mu)/dt  *  Fg,  where  Fg  is  the  buoyant  force  and  m,  is  the  mass  in 
the  vortex.  Since  m  «  R^  in  two  dimensions  and  Ru  is  constant,  the 
separation,  R,  increases  linearly  with  time  so  one  has  the  scaling  that 

R  «  t.  (10) 

Since  the  rise  velocity  of  the  vortex  pair  varies  as  r/R,  the  net  rise 
distance,  z,  increases  logarithmically  (in  2-D)  with  time.  Therefore, 
one  can  write  the  relation 


3 

For  three  dimensions  the  mass  varies  as  R  ,  the  momentum  equation  reduces 

2  i 

to  dR  /dt  equals  a  constant  after  torus  formation.  Then  R  <*  t*  and 

dz/dt  «  t‘*  and  unlike  Equation  10,  R  <=  z. 

The  two  dimensional  relations  (10)  and  (11)  were  tested  to  determine 
If  the  code  produced  the  correct  behavior  and  time  dependence.  Figure  1 
shows  the  result  of  the  test.  The  time  dependence  of  the  rise  rate  and 
vortex  widths  were  correct.  Figure  2  shows  the  time  evolution  of  the 
vortex  positions  as  the  bubble  evolves  into  a  torus.  The  evolution  was 
as  expected  and  showed  the  creation  of  the  vortex  pair.  Figure  3  displays 
the  corresponding  flow  fields,  note  that  the  vectors  are  renormalized  with 
Figure  3b  having  maximum  flow  ^  50  times  larger  than  3a. 

At  this  point  a  variety  of  experiments  were  performed  with  VINCE  to 
get  some  experience  with  the  various  nonlinear  interactions  one  might 
expect.  In  the  next  few  figures  a  series  of  time  dependent  snap  shots 
are  displayed  for  two  interesting  cases.  In  the  first  case  two  bubbles 
are  placed  on  a  vertical  line  as  shown  in  Figure  4a.  They  are  slightly 
more  than  one  radius  apart  at  t  *  0.  One  sees  that  evolution  has  only 
the  crudest  relation  to  the  1  bubble  situation.  The  bottom  bubble  does 
not  expand  horizontally  as  far  but  is  considerably  elongated.  At  a  time 
slightly  beyond  the  last  time  shown  the  lower  vortex  pair  actually  breaks 
off  from  the  material  rising  into  the  upper  torus.  One  then  actually  has 
four  discrete  regions,  the  two  vortexes  from  the  bottom  bubble,  the  center 
section  of  the  bottom  bubble  and  the  upper  bubble;  far  different  from 
the  expected  linear  picture.  In  Figure  5  one  sees  the  flow  fields 
corresponding  to  4a  and  4d.  Note  there  is  an  increase  of  -v  100  in  the 
maximum  flow  field.  Also  note  that  there  appears  only  two  discernable 
vortexes  in  the  flow  field.  The  remainder  of  the  large  flow  is  an  upward 
column. 

The  next  experiment  was  the  examination  of  three  cylindrical  bubbles 
rising  together.  Each  bubble  was  placed  approximately  one  radii  apart. 

As  can  be  seen  from  Figure  6,  the  subsequent  evolution  deviates  sub¬ 
stantially  from  the  one  bubble  case.  The  outer  halves  of  the  two  outside 
bubbles  evolve  in  what  appears  to  be  a  normal  fashion.  However,  the 


Figure  4.  The  two  bubble  case  is  plotted  here.  Note  the  considerable 

distortion  of  the  lower  bubble  as  the  two  interact.  The  bottom 
vortex  rings  have  essentially  separated  as  vortex  particles  of 
opposite  sign  are  almost  superimposed.  At  a  later  time  the 
two  bottom  regions  do  separate. 
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Figure  6.  The  vortex  positions  for  the  3  bubble  case.  Note  the  overlapping 
of  the  vortex  particles  in  plot  d.  This  indicates  a  separation  of 
this  region  from  the  rising  region  as  vortex  particles  have  opposite 
signs. 


center  region  does  not.  In  fact.by  the  last  time  step  shown, the  vortex 
regions  of  all  three  bubbles  have  not  only  broken  off  from  the  main  body 
of  the  buoyant  regions  but  they  have  failed  to  rise  substantially. 

Actually,  as  the  vortex  pair  expanded, the  Inside  pair  condensed  and  created 
pairs  with  a  flow  pattern  of  opposite  sizes  from  a  rising  vortex  pair. 

In  Figure  7  one  sees  this  reversal  of  flow  pattern  which  then  works  against 
the  buoyant  force  and  thus  retards  the  rise  rate. 

Because  of  the  expected  need  for  higher  dimensionality, another  portion 
of  the  simulation  effort  was  the  consideration  of  generalizing  the  vortex- 
in-cell  algorithm  to  three  dimensions.  After  some  analysis  and  experi¬ 
mentation  we  have  decided  that  such  a  generalization  is  possible  but  com¬ 
plicated.  The  main  thrust  of  course  Is  recognition  that  vortex  rings  may 
be  three  dimensional, but  vortex  particles  must  still  be  two  dimensional 
objects  that  are  generated  on  the  surface  discontinuity  with  a  directional 
vector  tangent  to  the  surface.  This  will  require  some  complicated 
Interpolation  onto  the  eularlan  grid  but  we  feel  it  can  be  done. 

The  reason  for  wanting  a  3-D  version  of  VINCE  Is  to  further  investi¬ 
gate  the  level  of  nonlinear  Interaction  in  cases  like  the  3  bubble  example 
shown  earlier.  Further, as  one  goes  higher  up  in  altitude  one  would  like 
to  be  able  to  consider  the  effects  of  an  anisotropic  pressure  tensor.  At 
altitudes  where  the  magnetic  field  may  play  a  noticeable  role, one  of  the 
first  effects  might  be  the  creation  of  an  anisotropic  pressure  tensor 
which  leads  to  terms  In  the  pressure  equation  called  magnetic  viscosity 

5 

terms.  To  date, the  effects  of  such  off-diagonal  terms  has  never  been 
explored  In  any  hydro  or  MHD  code  In  the  H.A.N.E.  community. 

To  summarize  our  progress  in  the  simulation  area.  The  algorithm 
apppears  to  have  considerable  potential.  In  the  coming  years  it  will  be 
expanded  Into  a  3-D  code.  It  will  be  used  to  test  such  models  as  the 
linear  superposlton  algorithm  now  In  place  In  C/LAMP.  This  algorithm  uses 
a  Hills  spherical  vortex  formation  to  superimpose  flow  fields  from  other 
bursts  onto  the  burst  region  being  considered.  The  code  will  also  be 
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SECTION  3 

CONVECTION  CONCEPTS 


The  other  major  task  during  the  contract  period  was  to  consider  the 
applicability  of  microstructure  techniques  to  the  mid  and  low  altitude 
regimes.  This  formalism  essentially  provides  a  method  of  identifying  fine 
scale  structure  Intensity  and  location  In  nuclear  burst  flow  simulations. 

The  microstructure  technique  was  first  applied  to  the  high  altitude 

7-9 

structure  problem  by  Dr.  Joseph  B.  Workman  and  Dr.  Frank  Chu  of  BRA 
and  has  since  been  successfully  implemented  by  BRA  in  the  SCENARIO  code.10 

Microstructure  theory  is  basically  the  application  of  stochastic 
process  techniques  to  flow  problems  that  are  typically  solved  numerically 
in  order  to  predict  the  location  of  structure.11  The  theory  has  been  used 
In  the  past  to  develop  algorithms  for  predictions  regarding  plasma  structure 
In  the  late-time  high  altitude  multiburst  environment.10  The  theory 
develops  the  algorithms  which  produce  estimates  of  structure  parameters  at 
the  modest  expense  of  additional  complexity  and  storage.  The  theory  itself 
Is  flexible.  The  methodology  of  Its  application  to  arbitrary  flow  problems 
Is  fairly  well  formulated.  The  success  of  the  resulting  algorithms,  how¬ 
ever,  Is  a  function  of  the  necessary  underlying  assumptions  regarding  the 
statistics  of  the  structure  and  is  not  guaranteed.  In  each  particular 
application  the  validity  of  the  statistical  assumptions  and  the  resulting 
algorithms  must  be  tested. 

To  date,  our  understanding  of  the  Information  necessary  to  use  the 
microstructure  approach  at  low  and  mid  altitude  is  not  complete.  However, 

some  observations  of  the  likely  phenomenon  of  importance  has  provided 

2 

direction.  In  the  high  altitude  algorithm,  ng  and  nfi  were  found  to  be 
statistically  conserved  quantities.  The  question  is  what  would  correspond 
to  such  quantities  in  lower  altitudes,  where  charge  density  does  not  pro¬ 
vide  the  structure  but  rather  radiating  material.  First, we  noted  that  at 
late  times  after  the  shocks  were  out  of  the  system,  the  flow  was  very  likely 
to  be  Incompressible  as  is  assumed  at  high  altitude.  Second, once  it  is 
observed  that  even  torus  formation  is  a  form  of  Rayleigh-Taylor  and  that 
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most  of  the  structuring  instabilities,  i.e.:  Rayleigh-Taylor,  Kelvin- 
Helmholtz,  Saffman-Taylor  and  such,  are  represented  by  the  creation  and 
convection  of  vorticlty  and  the  vorticlty  once  created  is  then  a  dynamical 
fluid  element  until  a  diffusion  mechanism  Is  provided,  one  might  reasonably 

consider  vorticlty  as  the  quantity  to  consider  when  the  role  of  n  or 

2  ® 
n‘  Is  explored.  This  assumption  must  be  tested  and  to  do  this  a  program 

to  test  them  was  devised.  Ue  plan  to  resort  to  the  vortex- in-cell  code 

as  a  test  devise. 

As  stated  earlier,  the  vortex-in-cell  code  formulation  focuses  on 
the  interface  between  the  hot  rising  fireball  gases  and  the  surrounding 
cold  air.  It  makes  the  assumption  that  the  flow  is  incompressible  and 
that  the  hot  and  cold  gases  do  not  diffuse  into  each  other.  The  interface 
between  the  two  gases  is  then  defined  and  followed  with  vortex  particles. 
One  can  assume  that  vortex  particles  move  and  vorticity  grows  in  a 
deterministic  manner  according  to  the  following  three  equations: 
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where  Is  the  vorticity  of  the  ith  particle;  v  is  the  fluid  velocity 
as  determined  through  Equation  (14),  the  potential  equation;  g  is  the 
acceleration  due  to  gravity;  ap  Is  the  change  in  density  across  the 
Interface;  and  az^  is  the  difference  in  altitude  over  the  region  covered 
by  the  ith  particle. 

Microstructure  theory  can  be  used  to  formulate  an  algorithm  which 
addresses  the  interface  structuring  issue  by  extending  the  particle  code 
technique.  The  basic  idea  of  the  extended  algorithm  is  that  the  position 
and  value  of  each  particle  vortex  Is  a  random  value  and  these  random 
values  are  characterized  through  their  mean  and  variance  which  the  algorithms 
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predict.  The  structuring  of  fireball  Interfaces  in  the  azimuthal  direction 
(which,  for  Instance,  is  apparent  at  later  times  near  the  bottom  of  ORANGE) 

Is  quantified  through  the  variance  of  the  interface  position. 

The  features  of  the  algorithms  which  have  been  established  are  that 
the  variance  in  vorticity  and  vortex  position  will  contribute  through  the 
potential  equation.  Equation  (14),  to  variances  in  the  flow  velocity. 

Variance  in  the  flow  velocity  will  lead  to  variance  in  the  vortex  position 
through  the  convection  equation.  Equation  (12).  Variance  in  the  vorticity 
of  each  particle  is  created  by  the  variance  in  the  value  of  az  which  is 
a  component  of  the  position  variance.  Features  that  remain  to  be  established 
are  the  detailed  nature  of  these  equations.  A  crucial  point  that  requires 
resolution  is  what  approximations  are  to  be  made  in  the  convection  step  to 
account  for  correlation  between  the  position  variate  of  a  particle  and  its 
velocity  variate.  An  assumption  of  no  correlation  could  lead  to  an  over¬ 
estimate  of  the  actual  structuring.  This  issue  and  similar  issues  will  be 
addressed  through  experiments  with  high  resolution  runs  of  the  VINCE  code 
and  with  further  analysis. 

To  sumnarlze  this  aspect  of  our  work.  We  believe  that  the  vorticity 
may  play  the  role  as  the  primary  structure  quantity,  the  existence  or 
variation  of  which  will  indicate  structure.  Considering  the  statistical 
properties  of  this  quantaity  in  much  the  same  manner  as  the  the  properties 
of  nfi  was  considered  in  the  high  altitude  may  provide  the  necessary 
information  for  a  structure  algorithm.  These  assumptions  will  be  further 
tested  and  numerical  examination  of  them  seems  possible. 


SECTION  4 


CONCLUSION 

Under  contract  no.  DNA001-84-C-0018  the  following  tasks  have  been  com¬ 
pleted.  A  vortex-in-cell  code,  VINCE,  has  been  constructed  in  two  dimensions. 
It  has  been  tested  and  several  multiburst  like  calculations  have  been 
run.  The  results  are  intriguing.  Their  correspondence  to  the  real  multi¬ 
burst  situation  has  yet  to  be  calibrated  but  will  be  against  full-fledged 
hydrosynamic  calculations.  As  configured,  the  code  can  do  buoyantly  rising 
bubbles  or  cylinders  since  these  are  2-D  and  ballistic  rising  bubbles  as 
well.  It  is  extremely  fast  with  the  calculations  shown  in  the  report, 
taking  less  than  a  minute  of  CRAY-1  time.  This  speed  will  allow  a  large 
variety  of  parameter  studies  into  the  various  structuring  mechanism. 

We  believe  we  have  found  a  way  to  apply  the  microstructure  concepts 
to  the  mid  and  low  altitude  problem.  Further  work  on  the  statistical 
properties  to  be  convected  will  have  to  be  done  and  the  VINCE  code  will 
play  a  role  in  evaluating  the  adequacy  of  any  future  model.  If  our  con¬ 
cepts  prove  valid,  structure  and  convection  models  for  the  mid  and  low 
altitude  should  be  possible. 

Finally,  we  investigated  the  possibility  of  extending  VINCE  into 
three  dimensions.  After  considering  various  problems  it  appears  that  it 
can  be  done.  We  will  be  making  this  extension  in  future  research.  This 
will  allow  us  to  do  much  more  realistic  studies  of  the  flow  fields  resulting 
from  multiburst  situations.  These  calculations  will  begin  from  either  a 
bubble  like  form  or  may  use  the  results  of  hydrodynamic  codes  as  starting 
conditions. 
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